Parquet approximation for the 4x4 Hubbard cluster 
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We present a numerical solution of the parquet approximation (PA), a conserving diagrammatic 
approach which is self-consistent at both the single-particle and the two-particle levels. The fully 
irreducible vertex is approximated by the bare interaction thus producing the simplest approximation 
that one can perform with the set of equations involved in the formalism. The method is applied 
to the Hubbard model on a half-filled 4x4 cluster. Results are compared to those obtained from 
Determinant Quantum Monte Carlo (DQMC), FLuctuation EXchange (FLEX), and self-consistent 
second-order approximation methods. This comparison shows a satisfactory agreement with DQMC 
and a significant improvement over the FLEX or the self-consistent second-order approximation. 

PACS numbers: 71.10.-w, 71.27.+a 



INTRODUCTION 



Over the past 50 years, many diflPerent techniques have 
been devised and employed to study strongly correlated 
electron systems. Unfortunately, advantages of the suc- 
cessful attempts were usually outweight by their limita- 
tions. Recently, because of the progress in computer tech- 
nology, complex diagramatic approaches have received 
increased attention. Although Baym and Kadanoff's $ 
derivability [l|, Q does not guarantee the physical valid- 
ity of a theory, their framework enables the generation 
of conserving approximations which are guaranteed to 
satisfy a variety of Ward identities. For these reasons, 
the FLuctuation EXchange (FLEX) approximation 0, 01 
has been intensively studied over the years. Its major 
disadvantage however is that it represents a conserving 
approximation at the single-particle level only. Thus, the 
physical validity of the approximation appears to be ques- 
tionable as the vertices are either overestimated or un- 
derestimated and the PauH exclusion principle is not re- 
spected properly 0|. In contrast, the parquet formalism 
[y| introduced by de Dominicis et al. in 1964 is a con- 
serving approximation which is self-consistent also at the 
two-particle level and one may hope that it resolves at 
least some of the limitations FLEX has. Unfortunately, it 
has extremely complicated structure and was, apart from 
applications to the Anderson impurity model and the 1- 
D Hubbard model with small system size 0, Q, hitherto 
also computationally out of reach. To circumvent this 
limitation. Bickers et al. introduced the so-called pseudo- 
parquet approximation j3| which attempts to improve on 
the FLEX without introducing the complexity of the full 
Parquet equations. But this approach fails to properly 
address the full frequency and momentum dependence of 
the scattering processes. Only very recently, due to the 
great advance of the parallel computing and the tremen- 



dous increase in computer memory, has it become possi- 
ble to fully solve this approximation for the first time. 

The paper is organised as follows. In part I we present 
the formalism and the resulting equations. In part II, we 
discuss the algorithm and the numerical difficulties that 
arising. In part III, we present first results obtained from 
the parquet approximation (PA) for the 2-dimensional 
Hubbard model and their comparison to other conserving 
approximaton methods such as FLEX and self-consistent 
second-order approximation (SC2nd). As a benchmark, 
we compare these results against the Determinant Quan- 
tum Monte Carlo (DQMC) which provides a numerically 
exact result. 



II. FORMALISM 
A. Vertex functions 

Standard perturbative expansions attempt to describe 
all the scattering processes as single- or two-particle 
Feynman diagrams. In the single-particle formalism 
the self-energy describes the many-body processes that 
renormalize the motion of a particle in the interacting 
background of all the other particles. In the two-particle 
context, with the aid of the parquet formalism, one is able 
to probe the interactions between particles in greater de- 
tail using the so-called vertex functions, which are matri- 
ces describing the two particle scattering processes. For 
example, the reducible two-particle vertex F^''(12;34) 
describes the amplitude of a particle-hole pair scattered 
from its initial state |3, 4) into the final state |1, 2). Here, 
i = 1,2,3,4 represents a set of indices which combines 
the momentum ki, the matsubara frequency itUm and, if 
needed, the spin ai and band index m^. 

In general, depending on how particles or holes are in- 
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FIG. 1: (color online) Different classes of diagrams; the solid 
line represents the single-particle Green's function and the 
wavy line represents the Coulomb interaction: here we use 
the p-h horizontal channel for illustration, (a) Reducible di- 
agrams: can be separated into two parts by cutting two hor- 
izontal Green's function lines, (b) Irreducible diagrams: can 
only be separated into two parts by cutting two Green's func- 
tion lines in the other two channels, (c) Fully irreducible dia- 
grams : can not be split in two parts by breaking two Green's 
function lines in any channel. 

volved in the scattering processes, one can define three 
different two-particle scattering channels. These are the 
particle-hole (p-h) horizontal channel, the p-h vertical 
channel and the particle-particle (p-p) channel. For the 
Hubbard model, the spin degree of freedom further di- 



vides the particle-particle channel into triplet and singlet 
channels while the particle-hole is divided into density 
and magnetic channels. 

One can further discriminate the vertices according to 
their topology. Starting from the reducible vertex F in- 
troduced above, we may define the irreducible vertex T 
corresponding to the subclass of diagrams in F that can 
not be separated into two parts by cutting two horizontal 
Green's function lines. Similarly, the fully irreducible ver- 
tex A corresponds to the subclass of diagrams in F that 
can not be split into two parts by cutting two Green's 
function lines in any channel. An illustration of these 
different types of vertices is provided in Fig. [TJ 

The Pauli exclusion principle produces the so-called 
crossing symmetries which in turn yield relationships be- 
tween these vertices in the different channels. This en- 
ables us to reduce the independent channels defined for 
the theory to the particle-particle and the particle-hole 
horizontal channels. 



B. Equations 

The parquet formalism assumes the complete knowl- 
edge of the fully irreducible vertices and provides a set of 
equations which are self-consistent at both the single- and 
two-particle levels. The connection between the single- 
and two-particle quantities is through the Schwinger- 
Dyson equation which connects the reducible vertex F 
to the self-energy E. It is an exact equation derived from 
the equation of motion and has the following form: 



S(P) - {G{P')G{P' + Q)G{P - Q){F4Q)p^Q,p, - Fr,,iQ)p^Q,p,) 

P'.Q 

+G{-P')G{P' + Q)G{-P + Q){Fs{Q)p-Q,p> + Ft{Q)p^Q.p,)} (1) 
I 



where G is the single-particle Green's function, which 
itself can be calculated from the self-energy using the 
Dyson's equation: 

G-\P)=G^\P) - i:{P) (2) 

Here, the indices P, P' and Q combine momentum k and 
Matsubara frequency iw„, i.e. P = (k, iw„). 

The reducible and the irreducible vertices in a given 
channel are related by the Bethe-Salpeter equation. It 
has the following form: 

Fr[Q)p,P' = Tr{Q)p^P' +<^r{Q)p^P' (3) 
Fr'{Q)p,P' = Tr.{Q)p,p. +-^r'{Q)p^P' (4) 

where r ~ d or m for the density and magnetic channels 
respectively and r' = s or t for the singlet and triplet 
channels, and we are using the vertex ladders which are 



defined as: 

MQ)rp' = Y.^-^'^^p^P"^o''iQ)p"^riQ)p",p'(5) 

P" 

^r'{Q)p.P' = ^Pr'(g)RP"Xo''(Q)P"rr'(g)p",p(6) 
P" 

Xo is the direct product of two single-particle Green's 
functions and is defined according to the particle-particle 
or the particle-hole channel. 

In a similar manner, the irreducible vertex and the 
fully irreducible vertex are related by the parquet equa- 
tion. This set of equations expresses the fact that the 
irreducible vertex in a given channel is still reducible in 
the other two channels. The parquet equation has the 
following form in the different channels: 
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TdiQ)pP' = AdiQ)pp' -^MP' -P)p.P+Q-l'^rniP' -P)P.P+Q 

+ i*,(P + P' + Q)^P^Q,^P + ^^t{P + P' + Q)-P-Q,-P (7) 

r™(Q)pP' = Am{Q)pP' - \^d{P' - P)P,P+Q + ^'^v^{P' - P)PP+Q 

- + P' + Q)-P-Q,-P + ^-^tiP + P' + Q)-p-Q.,-p (8) 

rs{Q)pP' = MQ)PP' + l'^diP' - P)-P'.P+Q - l^miP' - P)-P',P+Q 

+ i$rf(P + P' + Q)-p,,-P - \^m{P + P' + Q)-P'.-P (9) 

Tt{Q)pP' = At{Q)pP' + ^-i>d{P' - P)-P'.P+Q + l'i>rn{P' - P)-P'^P+Q 

- i$d(P + P' + Q)-P^^-p - i$™(P + P' + Q)-P',-P (10) 



The Bethe-Salpeter and parquet equations are also ex- 
act and derived from the categorization of the Feynman 
diagrams. 

The above discussion of the structure of the parquet 
formaHsm is far from being exhaustive and is merely in- 
tended to make the paper reasonably self-contained. For 
a more detailed description of the parquet formalism, we 
refer the reader to Bickers et al. Our actual goal is 

to numerically solve these equations self-consistently for 
the Hubbard model on a two dimensional cluster. The 
algorithm for this solution is described in the next sec- 
tion. 



III. ALGORITHM AND COMPUTATIONAL 
CHALLENGE 

The set of equations dicussed above are solved self- 
consistently as illustrated in the self-consistency loop in 
Fig. [21 One starts with a guess of the single-particle 
Green's function or self-energy. This can, for example, 
be taken from the second-order approximation. The re- 
ducible and the irreducible vertices are also initialized 
with the bare interaction. The self-consistency loop can 
then be described as follows: 

(i) first we calculte the bare susceptibility xo which is 
given by the product of two Green's functions 

(ii) next this bare susceptibility is used to calculate F 
through the Bethe-Salpeter equation 

(iii) we then proceed with updating the irreducible ver- 
tices F by solving the parquet equation. 
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FIG. 2: (color online) Schematic illustration for the differ- 
ent steps in solving the parquet approximation equations self- 
consistently. 

This step requires the input of the fully irreducible 
vertex A. In the context of the parquet approxima- 
tion which we study here it is taken to be the bare 
interaction. It however can also be extracted from 
some more sophisticated methods. 

(iv) it is followed by a calculation of the new F through 
the Bethe-Salpeter equation 

(v) this value of F is then used to update the self- 
energy through the Schwinger-Dyson equation 

(vi) the Dyson's equation is solved for the Green's 
function G. 

This loop is repeated until convergence of the self-energy 
E is achieved within a reasonable criterion. 

Unfortunately, this loop becomes unstable when the 
strength of the Coulomb interaction is increased or the 
temperature is lowered. As we beheve that this instabil- 
ity is purely numerical in origin and related to the itera- 
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tive nature of the algorithm, we have to extend the above 
scheme to account for this problem. For example, one 
possibility is to start with an overestimated self-energy 
and to damp it along with the irreducible vertex between 
two iterations according to: 



r 



, + (1 - ai)T,oid 
+ (1 - a2)Toid 



(11) 
(12) 



where ai and 012 are some damping parameters. 

Another possibility is to rewrite the coupled Bethe- 
Salpeter and parquet equations in the form /(x) = 
and apply a variant of a Newton's root searching method. 
Then we can take advantage of the existing linear solvers 
such as BiCGS GMRES [lH or the Broyden algo- 
rithm 

One major advantage that the parquet formalism has 
over Exact Diagonalization (ED) or Quantum Monte 
Carlo (QMC) is that it scales algebraically with the vol- 
ume of the system in space-time for any choice of pa- 
rameters including those that lead to a sign problem in 
QMC. The most time-consuming part of the formalism 
is the solution of the Bethe-Salpeter and the parquet 
equations, where the computational time scales as 0{nf) 
where rit = ric x nf, Uc being the number of sites on 
the cluster and n/ the number of Matsubara frequencies. 
Although the scaling is better than that of ED or QMC 
when the sign problem is severe, one can see that the 
complexity quickly grows beyond the capacity of usual 
desktop computers with incrasing system size, and large- 
scale supercomputer systems have to be employed. 

Our parallel scheme and our data distribution are 
based on the realization that the Bethe Salpeter equation 
is the most time-consuming part of our calculation. One 
can easily see that it decouples nicely with respect to the 
bosonic momentum-frequency index Q. This enables us 
to distribute the vertices across processors with respect 
to this third index and to solve the Bethe-Salpeter equa- 
tion with a local matrix inversion. However, this storage 
scheme puts a limit on the size of the problem that we 
can address. For a node with 2 GBytes of memory, the 
maximum value of nt that we can use if our variables are 
complex double precision is about 2500. 

Unlike the Bethe-Salpeter equation, one can readily 
observe that the parquet equation does not decouple in 
terms of the third index. Solving this equation requires 
a rearrangement of the matrix elements across proces- 
sors and this is the communication bottleneck in the al- 
gorithm. The rearrangement is necessary to obtain the 
form of the vertex ladder (f> or that is required in the 
parquet equation. For instance, in the d channel, we need 
<& (P — P')p pj^Q. This form of the vertex ladder is ob- 
tained by employing the three-step process described in 
the following equations: 



*(g)p,p, 
*(Q)p,p-p' 



^iP-P')p.Q 



(13) 
(14) 
(15) 
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FIG. 3: (color online) Single-particle Green function G (r) for 
the three diagrammatic approaches and the DQMC. For this 
temperature, the PA result (solid line) looks very close to the 
DQMC one (symbol solid line) as compared to SC second- 
order (dashed line) or FLEX (dash-dotted line). 



The first step in this transformation only moves data 
locally in memory. This does not require much time. The 
second step is actually just a 2D matrix transpose but 
with matrix elements spreading on many nodes. This 
is where communication across nodes is required. It is 
achieved by using the standard Message Passing Interface 
(MPI) collective directives The final step is also 

local and can equally be done very fast. 



IV. RESULTS 

In the following section, we will show the PA results for 
a 4 X 4 Hubbard cluster at half-filling. The calculations 
are done for U = 2t and different temperatures. The cal- 
culations are performed for a finite number of Matsubara 
frequencies [13|. However, for the observables we calcu- 
lated, such as the local moment and magnetic susceptibil- 
ity in Fig. m and Fig. [H we performed an extrapolation 
to an infinite number of frequencies so that the cutoff 
error in frequency is minimized. To see how good PA 
works for the lattice model, we use the DQMC result as 
the benchmark. In the DQMC calculation, Ar = 1/12 is 
used and the combined statistical and systematic errors 
are smaller than the symbols used. To further compare 
PA to other approximations, FLEX and self-consistent 
second-order results are also included. 



A. Single-particle Green function G (r) 

First, one can get a rough idea of how PA improves the 
accuracy of physical observables by comparing the single- 
particle Green's function from different levels of approx- 
imation. Shown in Fig. [3] are Gk (t) with k = (tt, 0) 
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FIG. 4: (color online) The inverse temperature dependence 
of local moment. Among the three diagrammatic approaches, 
the PA result comes closest to the DQMC one. 



calculated from the self-consistent second-order approx- 
imation, FLEX, PA and DQMC. The parquet result is 
significantly closer to the DQMC result than the second- 
order approximation and FLEX results as can readily 
seen from the figure. This confirms the intuition that one 
would get better results if the approximation is made on 
the vertex which is most irreducible. 



B. Unscreened local moment 

Next we present results for the local magnetic moment 
defined as 



(16) 
(17) 



where ha denotes the number operator for electrons of 
spin a. In the context of a conserving approximation, it 
can be re-expressed in terms of the self-energy and the 
single-particle Green's function as 



2T 



(18) 



where the trace sums over both the momentum and the 
frequency degrees of freedom. 

The results are shown in Fig. [H Among the three di- 
agrammatic approaches, the PA result comes closest to 
the DQMC one. If we look more carefully at the DQMC 
curve, we can find the existence of two humps. The hump 
at Ti ~ U/2, which is well reproduced by the PA, desig- 
nates the energy scale for the charge fiuctuation, and is 
directly related to the suppression of charge double oc- 
cupancy. The other hump beginning at r2 ^ ^ is related 
to the virtual exchange interaction, J, between nearby 
spins. It is believed to be related to the synergism be- 
tween the development of the long-range antiferromag- 
netic correlation and enhancement of the local moment. 
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FIG. 5: (color online) Uniform susceptibility calculated for 
different methods as a function of inverse temperature. While 
at the high temperature region, all the diagrammatic method 
results come close to the DQMC result, the PA shows its 
advantage clearly in the low temperature region. 



As a result, a pseudogap is opened which increases the 
entropy of the system [l^, The magnitude of T2 

can be estimated by noticing J = At'^/U for the strong 
coupling limit and t exp (— 27ri/J7) in the weak coupling 
limit [H, . Therefore it basically interpolates between 
these two limits for that U = 2t is in the intermediate 
coupling regime. This hump is not well captured by PA. 
The increasing importance of envelop-shape diagram con- 
tribution 0, H] not included in PA is responsible for this 
deviation in the low temperature region. 



C. Uniform susceptibility 

Finally, we look at the uniform magnetic susceptibility 
which is defined as 

Xrnag(0,0) = J^^T (fr S , (t) S , (0)) (19) 

= ^{S!) (20) 
with magnetic moment defined as 



SAr) 



N ^ 



("r,T (t) - nr,i (r)) 



(21) 



The Xmag from different approaches are presented in 
Fig. m The uniform magnetic susceptibility calculated 
from DQMC follows a nearly linear dependence on /?. 
This mimics closely the Curie- Weiss law of weakly inter- 
acting moments and implies that the dominant effect in 
the system is the short range magnetic fiuctuation. This 
is consistent with the f3 dependence of the local moment 
presented in Fig. [H As the temperature still dominates 
over the spin energy scale of the system, it suppresses the 
long range fiuctuation. 



6 



From this figure, the improvement of PA over the other 
two approximations is also easy to see. Similar to the 
local moment, the difference between results from PA and 
DQMC at the low temperature region can be explained 
by the omission of envelop-shape diagrams in PA. 

V. SUMMARY AND OUTLOOK 

We have presented the parquet formalism, PA method 
and in particular the implementation we use to solve 
large-sized problem. The preHminary application of PA 
on the 4x4 Hubbard cluster shows that it can yield bet- 
ter results than the self-consistent second-order or FLEX 
calculations. This is the first step in our work, next 
we are going to use the parquet formalism in the so- 
called Multi-Scale Many-Body (MSMB) approach 
Within MSMB, correlations at different length scales are 
treated with different methods. The short length scales 
are treated explicitly with QMC methods, intermediate 
length scales treated diagrammatically using fully irre- 
ducible vertices obtained from QMC and long length 
scales treated at the mean field level. Note that in this 
approach the fully irreducible vertex is approximated by 
a QMC calculation on a small cluster, while in PA it 
is approximated by the bare interaction. Therefor this 
approach should provide superior results to the PA. An- 



other advantage is that it can also avoid the exponential 
increase of the computational cost as the system size in- 
creases, and thus can take full advantage of the most up- 
to-date computer resources available. We will combine 
it with the Local Density Approximation (LDA) to gain 
some predictive power from the first principle electronic 
structure calculation. 
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